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(Dated: February 1, 2008) 

Abstract 

We present results for the static interquark potential, light meson and baryon masses, and light 
pseudoscalar meson decay constants obtained from simulations of domain wall QCD with one 
dynamical flavour approximating the s quark, and two degenerate dynamical flavours with input 
bare masses ranging from m s to m s /4 approximating the u and d quarks. We compare these 
quantities obtained using the Iwasaki and DBW2 improved gauge actions, and actions with larger 
rectangle coefficients, on 16 3 x 32 lattices. We seek parameter values at which both the chiral 
symmetry breaking residual mass due to the finite lattice extent in the fifth dimension and the 
Monte Carlo time history for topological charge are acceptable for this set of quark masses at 
lattice spacings above 0.1 fm. We find that the Iwasaki gauge action is best, demonstrating the 
feasibility of using QCDOC to generate ensembles which are good representations of the QCD path 
integral on lattices of up to 3 fm in spatial extent with lattice spacings in the range 0.09-0.13 fm. 
Despite large residual masses and a limited number of sea quark mass values with which to perform 
chiral extrapolations, our results for light hadronic physics scale and agree with experimental 
measurements within our statistical uncertainties. 

PACS numbers: ll.15.Ha, 11.30.Rd, 12.38.Aw, 12.38.-t 12.38.Gc 
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I. INTRODUCTION 



Quantum Chromodynamcs (QCD), a four-dimensional, asymptotically free gauge theory of 
interacting vector bosons and Dirac fermions, enjoys a unique position due to its mathe- 
matical consistency and, as the description of strong interactions, its relevance to current 
particle physics phenomena. However, due to the strength of the interactions in this theory, 
analytic calculations are not possible for many phenomenologically important quantities and 
numerical methods must be employed in the quest to achieve first principles results. Since 
Wilson's formulation of lattice gauge theory substantial progress has been made in these 
pursuits, although incorporating fermions on the lattice has been a particularly difficult en- 
terprise. Both Wilson's original lattice fermion formulation and the subsequent staggered 
fermion approach, at non-zero lattice spacing, break the full chiral symmetries of continuum, 
massless QCD. While these breakings vanish as the lattice spacing goes to zero, they lead 
to practical difficulties with operator mixing and a challenging extrapolation to the physical 
light quark masses. 

With the advent of domain wall fermions 2|, y, |4| and the related overlap fermions 5|, |6| , it 
has become possible to perform simulations at non-zero lattice spacing which possess good 
chiral and vector symmetries. When using domain wall fermions, at sufficiently small lattice 
spacings, the light-quark limit should behave as in the continuum, including the correct 
effects of topology, without any infrared pathologies. Chiral symmetry, which governs much 
of the physics of low-energy QCD, is mildly broken in a controllable way, yielding the correct 
number of light pseudoscalar mesons and protecting against the mixing of operators with 
different chirality. Also, as we will see, exact numerical algorithms, based on Monte Carlo 
techniques, can be employed with domain wall fermions for simulations with the three light 
quarks that are relevant to low-energy QCD. 

With domain wall fermions and exact algorithms for our simulations, we now have numerical 
methods which nicely complement the pristine mathematical pedigree of QCD. These meth- 
ods come at a cost. Current simulations with domain wall fermions require more computer 
power at a given lattice spacing than fermion formulations that do not fully realise chiral 
and/or flavor symmetry. This extra computational cost may be more than recouped by bet- 
ter control over the light-quark limit and potentially better scaling properties. Also, there 
are many calculations, particularly those involving nucleons and where operator mixing is 
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critical, such as those relevant for neutral kaon mixing, that may only be practical with 
fermion formulations that have good chiral properties. 

Simulations of QCD with two flavors of domain wall fermions were first performed some time 
ago on the QCDSP computer by the Columbia group g More recently, simulations at finer 
lattice spacings with better chiral symmetry were carried out by the RBC collaboration [8J. 
A first exploratory simulation of three-flavor QCD with domain wall fermions was also done 
by the RBC collaboration This work showed that dynamical domain wall simulations are 
practical, provided sufficient computer power is available. To increase the computer power 
available for QCD simulations, the RBC and UKQCD collaborations, in collaboration with 
IBM Research, the RIKEN Laboratory in Japan and the RIKEN-BNL Research Center, have 



10,01 



Both collaborations 



designed and built the QCDOC (QCD on a Chip) computers 
have a 12,288 node QCDOC, each with a sustained speed of 5 Tflop/s, making a serious 
investigation of QCD with 2+1 flavors of domain wall fermions possible. 
The primary objective of this paper is to establish a set of parameter values for simulations 
of 2+1 flavor QCD on QCDOC, using domain wall fermions and one choice from a class of 
improved gauge actions. The two light quark flavors (u and d) are given the same mass, mi, 
and we seek to make m; at least as low as one quarter of the strange quark mass, m s , which 
is held fixed close to its physical value. Clearly, we need to understand how our parameter 
choices affect the size of the residual chiral symmetry breaking in domain wall fermions. 
Since future simulations will require weaker couplings to test scaling, we also pay particular 
attention to the rate at which lattices decorrelate, using the evolution of topology as the 
most demanding measure. While small physical volumes will be used for the simulations 
presented here, once we have determined our parameters, volumes of 2.5 fm or larger, will 
be used. Since coarser lattice spacings make large physical volume simulations easier, we 
will also explore how domain wall fermions behave for coarser lattices. 
This paper is organized as follows. In Section [III we give our basic definitions and notation, 
the precise form of the actions used, and an overview of the simulations we have done. The 
work in this paper uses the exact Rational Hybrid Monte Carlo (RHMC) algorithm of Clark 

nn n 

and Kennedy [13|, H4J] and, for a few parameter choices, the inexact R algorithm [15J. We 
detail how we measure and extract masses, decay constants, the string tension and other 
observables in our simulations. In Section II III we describe our results for these quantities 
from simulations with a variety of gauge actions that differ in the gauge coupling, relative 
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size of the plaquette and rectangle terms, and light quark masses. Varying parameters in 
this space, keeping the physical lattice scale fixed, we measure and compare the residual 
chiral symmetry breaking and the evolution of topological charge. Since we envisage large 
simulations at several lattice spacings, we have also investigated the Iwasaki and DBW2 
gauge actions at weaker couplings. Our conclusions are given in Section IIV1 



II. CALCULATION DETAILS 



A. Lattice actions 



This paper reports on simulations with dynamical domain wall fermions for a class o 



proved gauge actions. Our notation, which we briefly review here, is the same as in |8j, [16 



lm- 



We denote points in four dimensional space-time by x and points in the fifth dimension of 
the domain wall formulation by s where < s < L s — 1 and L s is the extent of the fifth 
dimension. The partition function is given by 

/a 3 ft 3 

[du\ / ny*^] / Hi M 'pv/^pv^' 3 (i) 
J i=l i=l 

where the index i runs over the u, d and s quark flavors. The Pauli-Villars fields, <&pv,i are 
needed to cancel the bulk infinity that would be produced by the domain wall fermions as 
L s — > oo. In particular, the total action is 



S = S G (U) + S F {V, ^, U) + S PV (&, $, U). 



(2) 



The class of improved gauge actions we consider are of the form 



S G [U] 
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;i-8ci) £p[f/] 



X,fJ,U 



x;/i<u 



x;fi^u 



x,fj,u 



(3) 



where P[U] Xtfll> and R[U] x ,[w represent the real part of the trace of the path ordered product 
of link variables around the lxl plaquette and 1x2 rectangle, respectively, in the /x, v 
plane at the point x, and f3 = 6/g 2 with g the bare quark-gluon coupling. Different approx- 
imations to the renormalization group trajectory motivate two common choices for c\. 1] 



the Iwasaki action which sets c\ 
ci = -1.4069 



-0.331 



18 



19 



21 



22]. 



20| and 2) the DBW2 action which has 
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Jbr the fermion action in Eq. (jSJ), we use the domain wall fermion formulation of Shamir 



3fl, and Furman and Shamir 



4J. In our notation, the domain wall fermion operator D 



DWF 



for a fermion of mass rrif, is defined as 



D^,jM 5 ,m f ) = 5 S , S ,D"(M 5 ) + d,..,I)^,(m f ) 



(4) 



d 1 L(m 5 ) = - [(i - i,)u x j x+fl , x , + (i + -„)/;,/,. „. 

+ (M 5 -A)S xy 



(5) 



2 

mf 
2 



(1 - 7 5 )5 s+ i iS / + (1 + 75)5 s -i,y - 25 SjS / 

(1 - 75)^,L s -1^0,s' + (1 + 75)^,0^-1 



(6) 



D DWF is not Hermitian, but does satisfy 75-R5-DDWF75-R5 = -Ddwf' "where -R5 is a reflection 
operator in the fifth dimension. This, along with the transfer matrix formalism of 4J, suffices 
to show that det(.D DWF ) is positive for positive mass, allowing us rigorously to write the 
fermion action in a form suitable for simulations as 

3 



1/2 



(7) 



i=l 



where m« is the input bare quark mass for the ith light quark flavor. We only consider the 
case where all light quarks have the same value for the five- dimensional domain wall height, 
M5. The action for the Pauli-Villars fields is similar, except that the quark mass m/ is 
replaced by 1 to yield 

3 



S PV = J2®1 U4w(M>,1)Adwf(M 5 ,1) 



1/2 



(8) 



i=l 



.t should be noted that this is not the precise form of the Pauli-Villars action density given in 
4], but a variant introduced in j^j]. For the case of 2 dynamical flavors with the same mass, 
mi, integrating out the fermions and Pauli-Villars fields yields the following determinants 
to include in the generation of gauge fields 



det 



Dj, WF (M 5 , mi)D DWF (M 5 , mi] 



det 



d^ wf (m 5 ,i) j d D wf(m 5 ,i; 



(9) 
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This can be readily simulated by conventional Hybrid Monte Carlo by introducing pseudo- 
fermion fields for the numerator and using bosonic Pauli-Villars fields for the denominator, 
a 8 was done in Q. The m0r e recent 2-navo, Nations in fl U8 ed a S in g .e p 8 end„fen„io„ 
field to evaluate directly the ratio of determinants in Eq. This reduces the stochastic 
noise in the molecular dynamics evolution and speeds up the calculation, since a larger step 
size can be used. 

For 2+1-flavor simulations, the term in Eq. (Q has to be multiplied by 



(10) 



where m s is the strange quark mass. In this work, we have handled the fractional power in 
two ways. For most of the simulations we have used the exact RHMC algorithm {13I . \i\ . 



det 1 / 2 


-D D wf( M 5, m s )D DWF (M 5 , m s ) 


det 1 / 2 


Z^ WF (M 5 , 1)D DW f(M 5 ,1) 





For some of our initial simulations we used the R algorithm [15] , an inexact algorithm with 
finite step-size errors that easily handles fractional powers of the fermion determinant. At 
the time of this work, our implementation of the RHMC algorithm used separate fields as 
stochastic estimators for the numerator and denominator of Eq. ffTUl) . Code to stochastically 
estimate the ratio was being finished while these simulations were underway and is now in 
use. 

The RHMC algorithm allows us to simulate many decompositions of the same fermionic 
determinant, since det(M) = [det 1//n (M)] n . If we adopt the notation T){ra^) = 
D^ WF (M 5 ,m i )D- DWF (M 5 ,m i ) and, by convention, let every determinant appearing be rep- 
resented by a separate pseudofermion field, then our decomposition can be written as 

det^fPM] det 1 / 2 ^)] det^ 2 [V(m s )} 
det[P(l)] det 1/4 [£>(l)] det 1/4 [£>(l)] 

Thus, we used six pseudofermion fields: five associated with the RHMC, since five of the 
determinants in Eq. ffTTT) involve fractional powers, and one associated with conventional 
HMC. 



B. Details of the ensembles 

In our search for the optimal parameters for 2+1-flavor QCD with domain wall fermions, 
we have performed a large number of simulations for various values of (3 and c\ in Eq. ([3]) 
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with two different goals in mind: 1) to work at a fixed inverse lattice spacing (« 1.6 GeV) 
and see how the residual chiral symmetry breaking of domain wall fermions varies with c\, 
and 2) to then explore the weaker coupling behaviour at fixed C\. 

In earlier work using the QCDSP computer and the R algorithm, the RBC collaboration 
estimated that the inverse lattice spacing for 3-flavor QCD with domain wall fermions on a 
16 3 x 32 x 8 lattice with the DBW2 gauge action at (3 = 0.72 was 1.6-1.7 GeV Q. Since 
only a single dynamical quark mass was used and calculating the lattice spacing requires the 
light dynamical quark limit, they were only able to produce a rough estimate for the lattice 
spacing. We use the DBW2 action at (5 = 0.72 on a 16 3 x 32 x 8 lattice as our starting 
point. 



The RBC 



24( and CP-PACS 25|, [26] collaborations noted that the residual chiral symmetry 



;he Wilson 



171 ]. where 



breaking for domain wall fermions is reduced for the Iwasaki action compared to 
gauge action (c\ = 0) for quenched simulations. This was studied further in 
it was found that the DBW2 action markedly reduces residual chiral symmetry breaking 
on quenched lattices with a -1 ~ 2 GeV. For 2-flavor dynamical simulations with domain 



wall fermions, it was found in [24] that the Iwasaki gauge action was not much better than 
the Wilson gauge action for very coarse lattices with a -1 ~ 700 MeV. The recent 2-flavor 
simulations of QCD with domain wall fermions with a" 1 w 1,7 GeV Q] used the DBW2 
gauge action, but there was not sufficient computer power available then to test this choice. 
One of our goals in this paper is to pursue this question further to see how the gauge action 
choice affects residual chiral symmetry breaking, the tunnelling of topological charge and 
algorithmic performance. 

Thus, we explore three classes of gauge actions with different relative admixtures of the 
rectangle term in Eq. ([3]): Iwasaki, DBW2 and actions with even larger rectangle coefficients. 
All the ensembles were generated with a lattice size of 16 3 x 32 x 8. The Iwasaki and DBW2 
ensembles were generated with the RHMC algorithm and are described in Table [H All 
these ensembles have 2+1 flavors with the strange quark mass held fixed at approximately 
its physical value, and two or three values for the equal u and d quark masses, allowing 
for rudimentary extrapolations to the chiral limit. In addition to the quenched result that 
the DBW2 gauge action has excellent chiral properties 17], exploratory work [^tJ showed 
that a negative rectangle coefficient (c\ in Eq. (j3J)) has the effect of making the gauge 
fields smoother, thus rendering the residual mass smaller. We study this effect further, by 
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increasing the magnitude of ci, choosing (5 to keep the lattice scale fixed to approximately 
1.6 GeV. These ensembles were generated with 3 degenerate flavors using the inexact R 
algorithm, with the quark mass approximately that of the physical strange quark, and are 
described in Table [TTJ The three data sets at the bottom are at finer lattice spacings, which 
we generated to study the effect of this on the residual mass and the rate of change of 
topology. 

The RBC collaboration found that the combination of (3 and c\ that gives roughly the same 
lattice scale falls on the curve 

^ = -0.125 +A(a)P + B(a)p 2 , (12) 
Pp 

where (3r = c\ and j3p = 1 — 8ci. This is shown in Figure [H where early simulations with 
and 2 flavors of dynamical fermions and the simulations presented in this paper are included. 
The solid curve is a fit to the quenched simulations at an inverse lattice spacing of 2.0 GeV 
with different values of @r/ /3p. The ability to predict the lattice spacings in this manner 
allows us to do simulations at the desired coupling without the need for extensive searching 
in the parameter space. 

While the QCDOC computers provide some of the most powerful resources currently avail- 
able for lattice QCD, each parameter choice requires substantial computing resources. For 
example a 1,500 trajectory simulation takes a few weeks to generate on 1,024 nodes of QC- 
DOC, depending on the quark masses and algorithms used. In our study, the choice of 
parameters often depended on ensembles running at that time, so it was important to gen- 
erate ensembles as quickly as possible. We exploited naive parallelism and the availability 
of several smaller machines by spawning short Markov chains from the original chain, i.e. 
starting from a configuration in the original chain, a second distribution of random numbers 
(different from those in the original evolution) was generated and these were used to evolve 
a new branch. These branches were 'farmed' out to several machines in parallel. This had 
the advantage of increasing statistics to an acceptable level while reducing the 'wall clock' 
time. 

We used the decomposition in Eq. flTjJ for the (D, 0.72, 0.01/0.04) and (D, 0.72, 0.04/0.04) 
ensembles, whereas for all the other ensembles, including the (D, 0.72, 0.02/0.04) ensemble, 
we combined the two 1/4 power determinants in the denominator of Eq. ( TTTT) into a single 
1/2 power determinant and used one less pseudofermion field. 
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The parameters that enter into the RHMC algorithm control the accuracy of the rational 
approximation and the range of eigenvalues for which it is valid. Table UTTl gives the value for 
the parameters we used. The maximum and minimum eigenvalues, A max and A m i n , of T>{jrii) 
are used to determine the eigenvalue range for which the rational approximation has to be 
valid. The degree of the rational polynomial determines the accuracy of the approximation 
over this range. We use a more accurate rational approximation for the accept/reject step 
than for the molecular dynamics integration, since the accept/reject step removes any errors 
in the approximation arising during the integration. 

Table HVl gives values for (5H), (e~ SH ) and the acceptance for the ensembles generated with 
the RHMC algorithm. (e~ SH ) is equal to 1 within errors, indicating that the algorithm is 
working correctly. Also shown in Table [IV] are the ensemble averages for the plaquette, 
where typically the first 1000 trajectories of the ensemble were excluded from the averages. 



C. The static quark potential 

The static quark potential depends relatively weakly on the sea quark masses. Consque- 
quently, the chiral limit can be taken with reasonable confidence and provides an estimate 
of the lattice spacing that is relatively free of systematic uncertainty, compared to hadron 
masses. 

The static potential V(f) between a quark and anti-quark pair at relative spatial displace- 
ment r is calculated from the Wilson loop expectation value (W(r,t)), where 

(W(r, t)) = C{f)e~ v(f!)t + excited states. (13) 

Here (W(f,t)) is the average of the standard Wilson loop with spatial side of length f and 
temporal extent t. The static quark potential is then given by the ratio 

The time dependence in V(f, t) should disappear for sufficiently large t, so an important 
requirement is that V(r) be determined from a plateau seen in the quantity V(r,t) as t 
increases for fixed r. 

We compute lattice values for the Wilson loops (W(f,t)) using the method of Bolder et 



al. 



2n 



This approach evaluates all separations r, including those which do not lie along a 
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lattice axis. For such off-diagonal separations the Bresenham algorithm is used to determine 
the sequence of spatial gauge links that make up the spatial lines joining the two time 
segments of the Wilson loop. For both the on- and off-axis cases, these spatial links are 



constructed by APE smearing 



to improve the ratio of signal to noise. To be precise, at 



each step of APE smearing we replace each spatial link by 



U'(x) = P 



(U^U^x + v)Ul{x + /t)+ 



H = {1,2, 3}, (15) 



Ul(x - v)U^(x - i>)U v {x + /}-£)) 

where we include all spatial directions and P denotes a projection back onto SU(3). This 
process is then repeated iV smear 3( j times. The smearing coefficient and number of smearing 
steps are tuned to be (c smoari3 d, N smeaI >3( j) = (0.50,20 ~ 35), to maximise the overlap with 
the ground state of the Wilson loop, C(f). 

This approach was also used by the RBC collaboration ^_30] and an earlier description of 
some of the results presented here can be found in Ref. 3jJ. The results for V{f) obtained 
by this procedure were checked by an independent calculation using only on-axis loops and 
Chroma code 32|. While the errors were much larger for times t = 5 and 6, the results are 



consistent. 

Physical parameters are obtained by fitting the lattice value of V(r) to the function 



a 



V(r) = V --- + a\f\. 
\r\ 



Finally, the parameters a and a can be used to determine the Sommer scale 

r (mi) = 



33 



34] 



'1.65 -a) 



a 



(16) 



(17) 



for each gauge action and gauge coupling (3, in lattice units. The value of ro in the chiral 
limit can then be used to define the lattice spacing. To express this in physical units we 
take a = 0.5 fm/r . 



D. Topological Charge 

Continuum Yang-Mills gauge fields can be divided into classes characterized by an integer- 
valued winding number known as the Pontryagin index or topological charge: 

Q = ^ [ d 4 x e^rtr F^F^x) E Z (18) 
62ir z J 
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which is stable under smooth deformations of the gauge field. The Atiyah-Singer index 
theorem 35|, |36( predicts that the Dirac operator on a gauge background with topological 
charge Q will have at least Q exact zero modes in the zero quark-mass limit, 



Q = n + -n_, (19) 

where n + (n_) is the number of positive (negative) chirality zero modes. 
On the lattice the situation is much more complicated; gauge fields are far from smooth and 
the low-energy form of the Dirac operator can be distorted by the explicit breaking of the 
continuum symmetries. Nonetheless, when using domain wall fermions, zero modes of the 
Dirac operator are apparent and their numbers have been found to correspond well with 
the value of the topological charge calculated from the gauge field (by methods described 
below) fl6. I7I . 3?j]. Given the numerical cost in solving for the low-lying eigenvalues of the 
domain wall Dirac operator, we apply two gluonic definitions of topological charge: 

1. calculate the topological charge via Eq. ([TBI) using a classically 0(a 2 )-improved defi- 
nition of the field strength tensor_built from plaquette and rectangle clover-leaf terms 
(precise details can be found in 



2. the '5Li' definition of 



38 



39(, which combines the 1x1,2x2, 1x2, 1x3 and 3x3 
clover-leaf terms to give a classically 0(a 4 )-improved definition of the field strength 
tensor and, therefore, of the topological charge using Eq. ([18]) . 

Since both these methods expand about the classical continuum limit, they cannot be di- 
rectly applied to lattice gauge field configurations. Rather than deal with the subtle issue 
of constructing renormalised operators, we smooth the gauge field configurations so that 
expanding about the continuum limit is sensible. When using the first definition we follow 



401 ] in using APE smearing 



K(x) = p 



(1 - c„)u,{x) + -= {VMM* + fi)ut(* + P)+ 



Ut(x-0)U^x-0)U u (x + fi-u))] ,(20) 

with c smear = 0.45 in which the temporal link is also smeared in contrast to Eq. (1151) . We 
quote the value of Q after 20 steps, although we calculate up to 30 steps to check the 
stability of the extracted value. When using the second definition of topological charge, 
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we follow 



38 



391 ] and cool the configurations using the 5Li action. This takes the same 



combination of loops as used for the 5Li definition of the topological charge (although not 
clover- leaf symmetrized) to construct the action. Again, this definition is 0(a 4 )-improved 
about the classical continuum limit. It is also chosen so that the size of the instantons is 
invariant under cooling for instantons of size a few lattice spacings. We have used up to 50 
cooling steps. The results given in this paper are determined after 30 steps, beyond which 
the answers are stable. 

A third method we have used to monitor the topological charge is to calculate the fermionic 
operator (q^q) using a stochastic estimator. We only used a single stochastic source per 
configuration and, therefore, this measurement has large fluctuations when Q ^ 0, since the 
overlap of the noisy estimator with the small number of topological zero modes has a large 
variance. These fluctuations can be reduced by averaging measurements from nearby lattices 
in the Markov chain. At molecular dynamics time r = i, the smoothed value of (q^^q) is 
found by averaging over measurements in the range i — n/2 < r < i + n/2. We refer to this 
as a smoothing window of size n. Smoothing suppresses topological fluctuations that exist 
for only a short time in the evolution and should give a good estimator of the topological 
charge if the topology only changes significantly on a molecular dynamics timescale that is 
larger than n. The three methods give measurements of Q that agree extremely well. 



E. Hadron masses 



Hadronic correlation functions are constructed from quark propagators which take the fol- 
lowing form 

Sab{x,v) = (qA{x)q B (y)), (21) 

where the indices A and B represent different smearings of the quark field, q(x). All the 
smeared quark fields make use of the same construction of the four- dimensional quark field 
from the five- dimensional domain wall fermion field ty(x,s). In the case of the local (L) 
quark field, qii^), this takes the form 

q L (x) = P L *(x, a = 0) + P R V(x, s = (L s - 1)) (22) 
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where Pr/l — \ (1 ± 7s) are the chiral projectors. To construct the Coulomb gauge-fixed 
wall sources {W), we replace the quark field in Eq. fl2Tl) with the non-local field 



Qw 



{t)c = J2v(r,t)^q L (r,t) c , (23) 



where V(r, t) is the color matrix which transforms the spatial links at x = (r, t) into Coulomb 
gauge and c and c' are color indices. The third type of smearing (S) uses a hydrogenic 



wavefunction 41J. The quark field is convolved with a spatial smearing function after fixing 



the gauge fields to the Coulomb gauge as follows, 

qs(t) c = £ S (r, r ') V ( r '' *)c,c U{r\ t) d , (24) 
S(r,r')=e(- Ar /V8 r>ro , (25) 

where Ar is the modulus of the minimum distance between the center of the source, r , and 
r' (taking into account the periodic boundary conditions), and R is a tunable parameter. 
We give detailed comparisons of results for hadron masses obtained with different smearings 
in Section UTT1 

Correlation functions for mesonic operators are constructed as follows 

Cij(t,0) = ^e^(^(f,t)fi](0,0)) (26) 

X 

where 

n^t) =q(x,t)T i q(x,t) (27) 

and, for instance, 1^ = 75 for the pseudoscalar meson and Tj = 7^ for the vector meson. The 
spectrum can then be extracted by fitting these correlation functions at zero momentum to 
the form 

C(t) = A n {e~ rnnt + e- mn{T - t] ) (28) 

n 

where T is the size of the lattice in the time direction. For sufficiently large t, the ground 
state will dominate the correlation function. 

The standard baryon interpolating operator is composed of a local diquark operator and a 
spectator-like quark field: 

tti jk ,B( x ) = e abc [ql^{x)CTq h j(x)] q Cik (x), (29) 
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where T stands for one of the possible 16 Dirac matrices and C is the charge conjugation 
matrix. The superscript T denotes transpose and the indices abc and ijk label color and 
flavor, respectively. For the (J, J) = {\,\) baryons (the nucleon), an iso-scalar diquark is 
chosen, i.e {T = 1,75,757^}. However, only two of the three operators are independent, as 
they are related to each other through a Fierz transformation. In this work, the two nucleon 
operators are chosen to be 

Mbi = e a bc [ulC^d b ] u c , (30) 

^B 2 = tabc [UaCd b ] U c . (31) 

The intrinsic parity of these operators is defined by the parity transformation of the internal 
quark fields. The Qb 1 {&>b 2 ) operator transforms as a positive (negative) parity operator. 
However, this parity assignment is easily flipped by multiplication of the local baryon inter- 
polating operator by 75. Therefore, a two-point correlation function composed of either f^i 
or Qb 2 operators possesses both the positive- and negative-parity nucleon contributions. For 



details, see Ref. 
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431 ] . The operator is conventionally used in lattice QCD for the 



J p = l/2 + nucleon (N). However, it is also of interest to examine the J p = l/2~ nucleon 
(N*), so both Qb 1 and Qb 2 operators are utilized for (J, J) = (|, |) baryon spectroscopy. 
Taking the trace of the baryon two-point correlator with the relevant projection operator, 
P + = (1 + 74)/2, the two-point correlator in a finite box with (anti-) periodic boundary 
conditions is given by 

C B (t) = A v e- m ^ - A^e-™-^"^ , (32) 

where the parity of the forward propagating state 77 = + (— ) is the same as the intrinsic 
parity of the interpolating operator, ft Bl (^b 2 ), whilst the backward propagating state has 
the opposite parity. 



According to an intensive study of the nucleon excited states in Ref. |42j], the VLb 2 operator 
has a poor overlap with the nucleon ground state. Therefore, in this study, the masses of the 
positive and negative parity states were determined by a simultaneous fit to the following 
forms 

C Bl (t) = A + e~ m+t - A^e~ m - {T - l \ (33) 
C B2 (t) = A V m -'. 
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F. Residual mass 



While at short distances the domain wall fermion formulation is a five-dimensional lattice 
theory, at long distances and for large L s it is expected to appear identical to continuum 
QCD, with chiral symmetry broken only by the explicit mass term mj, introduced in Eq. (jH)). 
The deviations from this ideal behaviour can be easily described by a continuum, Symanzik 
effective Lagrangian. Because of the finite size of the fifth dimension, this effective La- 
grangian will contain explicit chiral symmetry breaking terms. While the coefficients of 
these terms are suppressed as L s increases, they are important because of their chiral prop- 
erties and because they are of lower dimension (three and five) than the dimension-six, 
chirality conserving, 0(a 2 ) terms: 

£ eff = L QC D(m f = 0) + (m f + m ICS )qq + c 5 qa^F^q. (34) 

The leading term is simply an additional mass term whose coefficient is called the residual 
mass, labelled m Ics . As is conventional, the normalization of m rcs is fixed by requiring the 
input parameter rrif and m ros to multiply the same qq mass operator as in Eq. (jMJ). In this 
paper we will ignore the effects of the Sheikholeslami-Wohlert, C5 term, which is suppressed 
by both large L s and one power of the lattice spacing. As will be seen below, methods for 
determining m rcs typically determine a quantity which depends on the input quark mass m/. 
This rrif dependence is a lattice artefact and represents the 0(a) or 0(a 2 ) ambiguity in our 
definition of m ies . A more careful treatment, which requires the analysis of the C5 term in 
Eq. (1341) and similar lattice artefacts in the quantities being used to compute m ies is beyond 
the scope of this paper. 

The residual mass can be computed from the additional contribution to the partially 



conserved axial current, Jg q , at the mid-point of the fifth dimension (J, [16|, |44|. As- 
suming that this mid-point contribution can also be described using the Symanzik ef- 
fective theory of Eq. ( 1341) . then to lowest order in the quark mass and lattice spacing 
J5 = m rcs J5 = m ros g7 5 t a g, where t a is a generator of the flavor symmetry. We can then 
compute m ros by averaging over time the quantity R(t), which is defined by 

R(t) _ (E* ./£,(*, *K(Q)) 

m ~ <£^(^M0)>- (35) 

Here 7r a (0) is a (possibly smeared) pseudoscalar interpolating field at t = 0. The minimum 
time used in the average need only be large enough to remove any contribution to the 

1(3 



correlators in R(t) from unphysical states, since m rcs should affect all physical states equally 
In the following we will refer to this average of R(t) as m' ics (mf) to explicitly display its 
dependence on m f and to allow a clear discussion of how we deal with this small, non-zero 
lattice spacing ambiguity in determining the constant m ros which appears in the effective 
Lagrangian of Eq. (I3"4"|) . 



G. Determining the pseudoscalar decay constant 

The pseudoscalar decay constant is defined by 

Z A (0\ q(x)Yl 5 q(x) |tt,p) = -ifp^e'^* (36) 

where Z A , the renormalization factor for the local axial vector current, can easily be deter- 
mined for domain wall fermions since Z A A a ^ = A a ^ and is the (partially) conserved axial 
current. In particular, we follow 16| and define 

C(t+l/2) = ^{AZ&t)*"^)) 

X 

Lit) = ^(AS(f,t)7r a (0,0)). (37) 

X 

We then calculate Z A from these correlators, using a combination that is free of 0(a) errors 
and minimizes 0(a 2 ) errors. This gives us the following explicit form for Z A . 

= 1 1 ( C(t + 1/2) + C(t- 1/2) 2 C(t + 1/2) \ 

A 1 + tmax - tmin J^. 2 \ 2 L ^ L ^ + + 1) J ' 1 ' 

v — ^min 

We have calculated fp in three ways. The simplest uses the local-local axial current in 
Eq. ([2S]), giving 

C £*>® = E<0|A)(*,*) 4>(0,0)|0> - A L A L oAo (e- mt + e~ m ^). (39) 

X 

The amplitude A A ^ A is determined from fitting the correlator to Eq. (139|) . which in turn 
using Eq. fl36l) yields fp from 

l2A LL 

f P = Z A J-^S.. (40) 
y nip 

In the second method, we use the axial Ward-Takahashi identity which, including the mid- 
point contribution as J£ q = m ICS Jg, gives 

m P Z A (0\A \7r) = 2(m f + m rcs )(0|P|7r) (41) 
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where P is the pseudoscalar density. From the local-local pseudoscalar correlator, we have 



Cg P {t) = ^2(0\P(x,t) P(0,0)|0> - A L P [ 



LL ( —mpt 



+ e 



-m P (T-t)- 



(42) 



which, combined with Eq. (1361) . gives 



fi 



2(rrif + m r 



'n ALL 



mp y mp 

The third method uses the local axial current pseudoscalar correlator 

(o|A)k>(7r|p|o> 



ciV*) = I>|4)(f,f) ^(o,o)|o) 



2m p 



e -mpt _|_ e -m P (T-t) 



)• 



(43) 



(44) 



The ratio of this correlator and the pseudoscalar density correlator is then fitted to the 
following form, 





7r)(7T 


P|0) 


(0|P 


7r)(7T 


P|0) 



to yield 



(0|A)k) 



tanh 



T 

mp(--t) 



(45) 



(46) 



(0|P|tt) 

Using the amplitude from the fit to the pseudoscalar density correlator, together with the 
value of Za, we obtain an expression for the pseudoscalar decay constant as 



fi 



2A 



p,p 



rap 



Aa ,p- 



(47) 



Comparing these different extractions allows us to probe the systematic error due to both 
mass and amplitude extraction, and the degree to which the chiral Ward identity is satisfied 
after shifting mf — > mj + m res . 



H. Autocorrelation length 



The integrated autocorrelation time is defined as 45j, |46| 



T , 



^ oo ^ oo 

^5>i(t) = ~+5>A(f) 



(48) 



where the autocovariance function, p^(t) defines the exponential autocorrelation time £ 



PA(t) 



PA(t) 



t^oo 



-t/Te: 



cxp 



(49) 
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and the autocovariance of an observable A is 



r A (t) = ((A s -(A))(A s+t -(A))). 



(50) 



The subscripts t and s label the Monte Carlo time and the outer average is over all pairs 
separated by t. It is standard to consider configurations separated by 2r A at to be statistically 
independent. 

In practice, we truncate the sum in Eq. (1481) at some finite value, t max , and define the 
cumulative autocorrelation time to be 



In the limit of sufficiently large t maxj this will be a good approximation to T A nt . In order to 
obtain reliable estimates for r exp and r^ um , autocorrelations should ideally be measured using 
ensembles containing many more configurations than the value of r A nt . Conversely, it must 
be noted that the measurement of the autocorrelation time is not sensitive to correlations 
over ranges of trajectories that are an appreciable fraction of the ensemble length, or greater. 

III. RESULTS 

Figures [2] and [3] shows example plots for the autocorrelation function and the integrated 
autocorrelation length for the plaquette and the pseudoscalar meson correlator at timeslice 
12. The statistical errors plotted for pA{t) and r cum were estimated using a jackknife proce- 
dure. In order to take into account the effects of autocorrelations in the error estimates for 
p{t) and r cum themselves, the original data for p(t) and r cum were grouped in bins of size b. 
Bin size was increased until the size of the jackknife errors stabilised to give the error bands 
shown on the plots. The final results can be found in Table IVl for both the DBW2 (3 = 0.72 
and Iwasaki — 2.13 RHMC datasets (the other ensembles are much shorter, and so the 
results are not given). 

A. Static quark potential 

The results for the values of the parameters Vo, a, u and r in Eq. ffl6l) and ffl7|) are given 
in Table IVljo^ ] . These are obtained by using the fitting range r G [V3, 8] and the choice 




(51) 



t=i 
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t = 5 in Eq. (|T4j) . The first error in tq given in the table is statistical and the second is an 
estimate of the systematic error in the fitting procedure. This estimate of systematic error 
is determined from the shift in the central value when the limits of the fitting range in r are 
swept through r min £ [\/2, a/6] and r max £ [7,9] and t is changed from 5 to 6. In Figure @] 
we show some example plateaux found for tq as t is varied between 3 and 6, suggesting a 
reasonably stable result and good suppression of excited states by t = 5. 
In Table IVHI we list the result for r that we obtain from a simple linear extrapolation to 
vanishing light sea quark mass for the four ensembles where we have two or more values for 
the light sea quark mass. Figure [5] compares two extrapolations for the DBW2 (3 = 0.72 
ensembles where we use either the two lightest mass values, mi = 0.01 and 0.02 or all three 
mi = 0.01, 0.02, and 0.04. 

B. Topological Charge 

The time histories of the plaquette (qq), (<?75<?) and the topological charge, are plotted 
in Figure [6] for the (D, 0.72, 0.01/0.04) ensemble. Measured values are plotted every 10 
trajectories (5 molecular dynamics time units). This ensemble has the smallest values for 
mi, and these time histories make our choice of 500 time units for thermalization appear 
reasonable. 

The topological charge in the bottom panel of Figure[6]was measured using gluonic method 1, 
as described in Section HlDI We can compare this with a fermion based definition of topology 
using ((775(7). To make this comparison, Figure [7] gives the evolution of a smoothed version of 
(qi§q) f° r the (D, 0.72, 0.01/0.04) ensemble, with smoothing windows of size 25, 50, 100 and 
200 time units. Comparing these evolutions with the unsmoothed ((775(7) evolution in the 
third panel of Figure [6] reveals noise is substantially reduced by smoothing. The relevance of 
the resulting signal can be seen by comparing with the evolution of the topology as measured 
from the gauge field. Figure [8] directly compares this estimate of the evolution of topological 
charge with the evolution of ((775(7), which has been smoothed with a window size of 50. 
The data are strongly correlated and show good agreement between topology as determined 
from smeared gauge links and from the Dirac operator. Since smoothing the data changes 
its normalization, we have rescaled the smoothed data such that its largest value is equal to 
that of the largest value of topological charge from the gauge field. 
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A comparison of 5Li and classically improved methods has been performed on the (I, 2.13, 
0.04/0.04) ensemble. Figure [9] shows the topological charge history from the two methods 
As can be seen they track each other quite well, although they disagree in places by as much 
as ~ 3 units of topological charge. It can also be noted that the 5Li method give results 
which are much better aligned with integers than the classically improved method. 
As a whole, while the topological charge is certainly not completely decorrelated from con- 
figuration to configuration, this set of ensembles are sampling the topological sectors quite 
well. This can also be seen in the histograms shown in Figure [TOl 

Figure dH shows the gluonic measure of the topological charge for the large rectangle simula- 
tions at fixed lattice spacing. Although the topological charge is sampled well between —15 
and 15 for the (D, 0.72, 0.04/0.04, R) ensemble (top panel), as we make c\ more negative, 
the range of the topological charge we are able to sample tends to narrow. For the (C7.47, 
0.16, 0.04/0.04, R) ensemble (bottom panel), the topological charge fluctuates only in the 
range of ±5, and it takes a large number of trajectories to move between different topology 
sectors. The situation is worse at weaker couplings, where the topological charge evolution 
is much slower. This can be seen in Figure H21 where the top three panels have a rectangle 
coefficient of c\ = —3.57 and the bottom one has c\ = —2.3. In the absence of a dramatic 
improvement in the chiral properties of domain wall fermions, this is a compelling reason 
to avoid larger rectangle coefficients. As should be expected, these trends with rectangle 
coefficient and coupling also apply to the DBW2 and Iwasaki actions. Figure [13] shows - 
from top to bottom - representative topological charge histories for the DBW2 action with 
P = 0.764 and (3 = 0.78, and the Iwasaki action with (3 = 2.13 and (3 = 2.2. 

C. Hadron masses 

Hadronic correlation functions were constructed from quark propagators which were smeared 
at the source and local at the sink, as described in Section III El We use the following 
notation to describe the combinations of smearing: SiFi — S2F2, where Sj denotes source 
and Fj sink for the quark propagators % = {1,2} in the meson correlator. In particular, four 
combinations were computed: LL-LL, SL-LL, SL-SL, and WL-WL. The first combination 
has the attractive property that, as both meson interpolating operators are the same, the 
ground state and all the excited state contributions to the correlation function are positive. 
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The effective mass, 



M eS (t) = log 



C(t) 



(52) 



C(t+l)_ 

which is used to determine when the ground state is dominating the correlation function, 
then approaches a plateau from above and can be used to determine unambiguously a fitting 
range for a single exponential function. Smeared and wall sources have the advantage that 
the plateau, and hence the fitting range that can be used, starts earlier. 
Figure shows the pseudoscalar effective mass for four different source/sink combinations 
as determined from the pseudoscalar-pseudoscalar correlator for the (D, 0.72, 0.01/0.04) 
ensemble, while Figure [15] shows a similar plot for the vector meson effective mass. In both 
cases, the mesons are constructed from degenerate mass quarks with m va i = 0.01. These 
plots use a radius of 3.0 for the smearing function (R in Eq.( 1231) ). which was found to be 
approximately optimal for both the pseudoscalar meson and the vector meson. This can 
be seen from Figure [161 which shows a comparison of the vector meson effective mass for 
various different smearing radii on the (I, 2.2, 0.02/0.04) ensemble with a valence quark 
mass of 0.04. When multiple sources are available, we have also extracted the mass by 
simultaneously fitting a pair of correlators with different sources to the ground and first 
excited state. The advantage of this approach is that the systematic uncertainty in the 
ground state mass arising from the choice of fit range is reduced. In summary, with careful 
choices of fitting ranges all the sources used in this work give consistent results. 
The extracted pseudoscalar meson masses for the cases where valence and sea quark masses 
are equal, i.e. unitary data, are given in Table lYTiTl for all of the RHMC ensembles. For 
the DBW2 gauge action at j3 = 0.72 we also have data for which the valence quarks are not 
degenerate with the sea quarks. These results are given in Table IIX1 The corresponding 
results for the vector meson masses are given in Tables |X] and IXII respectively. 
As with the construction of the meson operators, smeared quark fields are used to construct 
improved nucleon interpolating fields. In particular, all the data presented were extracted 
using correlators in which one interpolating operator is constructed from local quark fields, 
while the other interpolating operator is either constructed from smeared quark fields (de- 
noted SL-SL-SL), or gauge-fixed wall sources (denoted WL-WL-WL). Figure [TTl shows a 
typical effective mass plot for the nucleon. The backwards propagating, negative parity 
state in the standard nucleon correlator has been reflected about the middle of the time 
axis. Both states are shown together in Figure [T71 where solid symbols represent data ob- 
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tained from the WL-WL-WL correlators while open symbols correspond to the SL-SL-SL 
correlators. Fits to either type of smearing produce the same nucleon mass within the 
measured uncertainties. 

The mass of the negative-parity nucleon is less well determined due to the poor signal and 
has some dependence on the type of smearing chosen. Moreover, the physical volume of these 
lattices is likely to be too small to extract accurately orbitally excited baryon states such as 
the N*. However, the data presented can be used in combination with larger volume runs 
to estimate the size of the finite volume effects on the remaining spectrum. The extracted 
values of the nucleon mass are collected in Table Kill The three fit ranges given are for the 
forwards propagating state (iV + ), the backwards propagating state (A r _), and the negative 
parity correlator (N'_). 

D. The residual mass 

The ratio R(t) given in Eq. (135 j) is shown in Figure fT8l for the ensemble (D, 0.72 0.01/0.04). 
Since the correlator is time symmetric, the data shown is the average of the two halves of 
the lattice. This ratio should be constant as long as the separation in time is large enough 
that the effects of non-physical states/non-locality are eliminated. For the LL-LL correlator 
this is from t > 7 in the figure, but it is earlier for the WL-WL correlator as the overlap of 
the smeared operator with the physical states is larger in this case. As discussed above in 
Section III F\ we refer to this average of R(t) as m' es (mj), to distinguish it from the constant 
m ros which appears in the Lagrangian of Eq. fl34|) . The values of rn' TBS { m f) fr° m the unitary 
data are given in Table KHII and from the non-degenerate data in Table IXIV1 
Table IXVl contains the results for m' res (mf) and the vector meson mass as the magnitude 
of the rectangle coefficient c\ is increased. In this case m' ics (mf) has been estimated by 
averaging R(t) from t = 4 to t — 16. The decrease in m rcs as |ci| increases is expected, 
since the non-locality introduced by the large rectangle term suppresses the fluctuations of 
the gauge fields on the scale of two lattice spacings, and the number of low-lying localized 
m o d e S is m ode rat e ly .educed 0. If a SmaU residual mass were the only consideration for 
our parameter choices, our goal could be accomplished by increasing the magnitude of the 
rectangle coefficient. However, we would forfeit our ability to tunnel between different topo- 
logical sectors, as shown in the previous section. It is for this reason that we concentrated 
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on the Iwasaki and DBW2 actions in our more extensive action studies. 



E. Za and the pseudoscalar decay constant 

The value of Za is determined from Eq. (138]) . using either LL-LL or SL-SL correlators. In 
the case when smearing is used, this corresponds to using a smeared-smeared pseudoscalar 
interpolating field in Eq. f[3"T|) . Since the axial current operator is unchanged, this still 
corresponds to an extraction of the renormalization factor for the local axial current, A ^. 
Figure dn] shows a comparison of two such extractions for the (D, 0.72, 0.02/0.04) ensemble, 
with good agreement being found. The values of Za for the unitary data are given in 
Table IXVll For the (D, 0.72, 0.01/0.04) ensemble we also have some partially quenched 
data, extracted using LL-LL correlators. The values are given in Table IXVIII 
The value of the pseudoscalar decay constant has been determined in three ways according 
to Eq. (l40l) . ( 1431) and ( 1471) . These results are presented for both the unitary and non- unitary 
data in Tables IXVIIIIIXXll Figure [201 shows an example of the third method (Eq. (T4T1) ) for 
the DBW2 (3 = 0.72 ensembles. 



F. Chiral extrapolations and results for phenomenological quantities 

To estimate values of phenomenologically relevant quantities requires extrapolating our data 
to the physical u and d quark masses. For the DBW2 j3 = 0.72 ensembles we have three sea- 
quark mass values, as well as non-degenerate valence quark masses, enabling us to explore 
the extent to which our data agrees with the predictions of leading order (LO) or next-to- 
leading order (NLO) (partially quenched) chiral perturbation theory (PQxPT). Preliminary 



work on this data can be found in 
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For the remaining DBW2 and Iwasaki ensembles, 
we only have two sea-quark mass values and so only linear extrapolations are possible. For 
the larger rectangle ensembles, only one sea-quark mass was generated, so valence only 
extrapolations are performed. 

For each of these extrapolations we must include the effects of the residual quark mass 
which should be treated as an addition to the explicit quark mass rrif. Thus, the chiral 
limit should be defined as that value of mj for which m/ + m rcs = 0. Our values for the 
residual mass, m' es (m/), determined from the midpoint term in the axial current divergence, 
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themselves depend on mj, a lattice artefact. We treat this nif dependence in two different 
ways where we have three sea-quark masses. Firstly we extrapolate linearly in mj to define 
m rcs = m' res (mf = 0) for the unitary data, giving m ies = 0.0106(1). This is shown in 
Figure [2T] Secondly we construct the quark mass 

m q = m f + m' iea (m f ) (53) 

and extrapolate hadronic quantities versus m q to m q = 0. For data where we have two sea- 
quark masses we follow the second procedure. These two procedures adopt slightly different 
choices for the residual mass: m ros = m' cs (mj = 0) and m rcs = m' cs (mj = — m ies ). The 
difference between these two approaches is at most a few percent and reflects a non-zero 
lattice spacing error that we do not attempt to control in this paper. 

Concentrating on the DBW2 (3 = 0.72 ensembles, Figure [22] plots the square of the pseu- 
doscalar meson mass from the unitary data versus rrif. Lowest order chiral perturbation 
theory would suggest that this should be linear in the mass, and extrapolate to zero in 
the rrif = —m rcs limit. This is not exactly obeyed, mp extrapolates to —0.0144(24) at 
rrif = — m rcs . This is due to the relatively small size of the fifth dimension. Conversely, 
vanishing pseudoscalar meson mass occurs at a value of rrif = —0.0078(5). Whilst this 
represents a relatively large discrepancy with the definitions of zero quark mass used above, 
and is much bigger than the discrepancy in the definitions of m ros , it is itself a small effect 
when compared to the statistical uncertainties of other hadron masses and matrix elements 
extrapolated to the chiral limit. Ultimately, the uncertainty in the value of m ies is not the 
dominant source of uncertainty in the value of any results presented here. 
At NLO in PQ^PT the quark mass dependence of the pseudoscalar meson mass and decay 
constant has both a linear and non-linear component. For the lightest two DBW2 (3 = 
0.72 ensembles we have generated non-degenerate data for these two quantities, detailed 
in sections IIII CI and IIIIEI Consequently, we can attempt fits to NLO partially quenched 



chiral formula to this data. A preliminary study of these fits was first reported in [48j where 
nip and fp were fitted to the chiral formula independently. Here we update the results 
by performing combined fits to both nip and fp. While we observed consistency between 
QxPT and our data with valence quark masses as heavy as 0.03 in the independent fits of 



48j , this is not the case for the combined fits, shown in Figure [22] and [24] The fitted curves 



miss the data points badly at large quark masses. This is likely due to the fact that the quark 
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masses are so heavy that the next-to-next-to- leading order contribution becomes important. 
Given the limited statistics, large residual mass, relatively small volume and coarse lattice 



spacing, our study on the chira 
lattice artifacts are in progress 



fits is far from conclusive. Further investigations on these 
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For all the RHMC datasets, we extrapolate linearly in quark mass {m q ), the unitary data 
for the quantity mx as follows 

fn x (m q ) = A x + B x m q , (54) 

where mx = {mp, my, m^, m^*, fp}- Examples of these extrapolations for various datasets 
are shown in figures [2511281 

To predict physical quantities our strategy is as follows: The u and d quark masses are set 
to zero, and the lattice spacing is determined, either from r or by setting the vector meson 
mass in this limit to be equal to that of the physical rho meson. The s quark mass is then 
set from the physical kaon mass in lattice units and the pseudoscalar meson mass made 
from two quarks each with a mass of half the strange quark mass. Predictions for the values 
of the following quantities can then be made: {m^*, m^, m^*, fx}- Preliminary results 



were first presented in 50|, [51 ] . 



The results for the value of the vector meson in the chiral limit are listed in Table IXXIII 
The lattice spacing can then be determined from the physical rho meson mass. Whilst not 
ideal, as the rho meson can decay in nature and has a large width, and because we should 
extrapolate our data to the u and d quark masses, not the chiral limit, this method should 
be accurate enough for our purposes. The values of the lattice spacing determined in this 
way are also displayed in Table KXHl One can compare these values for the lattice spacings 
with those obtained from the potential by setting tq = 0.5 fm, in Table rVHl There is crude 
agreement, but the lattice spacing from the potential is systematically finer. Alternatively, 
we can estimate the value of r$ from the potential and the lattice spacing determined from 
vector meson mass. This is also displayed in Table IXXIII where the errors quoted have been 
obtained by adding the errors for r^/a and a -1 from the vector in quadrature. A crude 
estimate for the value of r in the continuum can be made by an average, weighted by the 
square of the errors: 

r = 0.554(21) fm. (55) 
This is equivalent to a constant continuum extrapolation and ignores the different systematic 
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uncertainties for each ensemble, as well as the issues described above for extrapolation of 
the vector meson mass, and so should be only considered as a qualitative, rather than 
quantitative result. The determinations of the lattice spacing from the potential and the 
vector meson mass are consistent, given the uncertainty in the physical value of r . In 
general, despite its problematic nature, we used the vector meson mass to set the scale, 
where physical units were required and ro for scaling analyses as the value of r^/a is properly 
controlled for each ensemble. 

For the larger rectangle data, where only valence chiral extrapolations are possible, we set 
the lattice spacing from the vector meson mass. The results are displayed in Table IXXIIII 
Extrapolating my with respect to the dynamical mass instead of the valence mass will 
decrease the value of am p by 10% to 15%, and hence increase the value of a' 1 by the same 
percentage. However, this is likely to effect all four datasets by a similar amount, so the 
valence chiral my is good enough to compare the relative lattice spacings. We see that they 
are approximately matched, making the comparison of the residual mass meaningful. 
The approximate value for the bare s quark mass is extracted by requiring that a pseu- 
doscalar meson made up of two quarks of mass m s /2 has the experimental kaon mass. In 
doing this, we are neglecting both next-to-leading order chiral perturbation theory effects 
and the masses of the u and d quarks. Table IXXIVI collects together the intercepts for the 
linear extrapolations, plus the extracted strange quark masses. These chiral extrapolations 
were performed on the multiple correlator fits tabulated in Table IVIIII As can be seen, the 
canonical heavy dynamical quark mass we have been using (0.04) is very close to the value 
of the s quark mass extracted from this method. 

The values for K*, the nucleon and its negative parity partner are shown in Table IXXVI 
The determination of ma* can be seen for the DBW2 ensembles in Figure [26j We examined 
the scaling behaviour of dimensionless ratios versus (a/ro) 2 . For the vector masses, this ratio 
was / vn p and is shown in Figure [29j Given the change in lattice volumes, the relatively 
large values of residual mass, and only two sea-quark masses, the scaling behaviour appears 
reasonable. 
The J parameter 



52J is defined as 

dmy 
dm 



J = m v^ZT ( 56 ) 



p 
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and is determined at the experimental ratio 

M K * 



1.8. (57) 



M K 

The linear fitting forms for the chiral extrapolations of the vector {my) and pseudoscalar 
(mp) masses implying that a plot of my versus mp will be a straight line for varying quark 
masses. Here we are working along the 'unitary trajectory' where m sea = m va \. Figure [30] 
shows a plot of my versus mp for the DBW2 (3=0.72 and Iwasaki (3=2.13 cases and, indeed, 
for the (3 =0.72 case, where there are three points, approximate linear behaviour is seen. 
The intersection of this line with my = 1.8ms (starred points) determines the reference 
value my which is to be multiplied by the slope to yield J. Figure [31] shows the value of the 
J parameter on all the datasets. Within the large errors, agreement with the experimental 
value is observed. 

Figure [321 shows the dependence of the baryon spectrum at the chiral limit, in dimensionless 
units, on the lattice spacing. As the size of lattice is fixed to be 16 3 x 32, the physical extent 
of the box in fm is smaller for the smaller (a/r ) 2 in Figure [32] The scaling behaviour of 
the spectrum appears promising despite the change in volume. A continuum extrapolation 
is not attempted. The N* mass at the finest lattice spacings is approximately two sigma 
smaller than that at the coarsest lattice spacing, which is consistent with the theoretical 
expectation that the N* becomes degenerate with the N in a small enough box 42[. This 
would also suggest that finite size effects may be beginning to affect the N for the finest 
ensembles. These finite size effects would tend to increase the mass of the iV 53j . With the 



current statistical resolution it is not possible to judge whether this is really the case. 



The Edinburgh plot 54J is shown in Figure [33] This is a useful way of comparing results 
for different actions without the need for any extrapolations of the data. Shown on the 
graph are the experimental ratios and the values obtained in the static quark limit, where 



the hadron mass is equal to the sum of the valence quark masses. 



he curve obtained from 



the phenomenological model for the hadron masses described in |55J is shown as a guide for 
the eye. It is remarkable that, even at relatively coarse lattice spacing, with a small fifth 
dimension and moderate chiral symmetry breaking, the data follows the phenomenological 
curve, albeit with relatively large statistical error. This very promising result suggests that 
a reliable chiral extrapolation could be performed with more data. 

The values of the pseudoscalar meson decay constant in the chiral limit for all datasets are 
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displayed in Table IXXVI1 Setting this equal to the physical value of f n provides another 
estimate of the lattice spacing. This is also tabulated in Table IXXVII The different methods 
for extracting the decay constant give consistent answers. For the Iwasaki /3 = 2.13 data set, 
the AA correlator gives slightly higher values than the other methods. However, at weaker 
coupling all the definitions agree, perhaps suggesting that this effect is a lattice artefact. 
Figure [341 shows the scaling behaviour of the dimensionless ratio Jk/ from the PP method. 
This shows excellent scaling, albeit with large errors and, moreover, is in agreement with 
the experimental value. 

IV. CONCLUSIONS 

Our criteria to be able to simulate 2+1-flavor QCD on QCDOC using domain wall fermions 
are as follows. The topological charge autocorrelation length should not be significantly 
greater than one hundred HMC trajectories. Too infrequent changes in topological charge, 
particularly on the small lattices we are using, would signal poor ergodicity. On the other 
hand, the greater the presence of dislocations which drive topological change, the larger the 
density of localised near- zero modes which contribute to m res . To reach the chiral regime 
under theoretical control m rea has to be small compared to the explicit u and d quark masses. 
Thus, a balance has to be struck between increasing the rate of local topological fluctuations 
and decreasing m res , i.e. the level of chiral symmetry breaking. Finally, to control both 
lattice artefacts and finite volume effects for light hadron physics, we need to simulate at 
lattice spacings in the range 0.09-0.13 fm and spatial extent of at least 2.5-3 fm. 
We find the Iwasaki action satifies all these criteria. While the DBW2 gauge action produces 
a smaller m ICS at a given lattice spacing, we find the sampling of topological charge, especially 
at the finer lattice spacing to be too slow. Increasing the coefficient of the rectangle further 
produces only a modest decrease in m rcs , while further supressing topology change. 
We have presented results for the the static interquark potential, light meson and baryon 
masses, and light pseudoscalar decay contants. We find that, even on our small volumes and 
small extent in the 5th dimension, with rather heavy sea-quark masses and crude approach 
to chiral extrapolation, both the Iwasaki and DBW2 actions, reproduce experimental values, 
albeit with large statistical uncertainties. The scaling behavior appears promising over the 
range of lattice spacings from 0.09-0.13 fm, except for the baryons, where it appears to be 
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spoiled by finite size effects for the smallest volume. 

We conclude that, using the Iwasaki gauge action, there is a range of values for the quark 
masses and lattice spacings such that the chiral regime for 2+1 flavor QCD with domain 
wall fermions on QCDOC is accessible, while maintaining control of lattice artefacts and 
finite volume effects for a range of light hadron physics. 
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TABLE I: 2+1-flavor ensembles and the mnemonics used to describe them. A * denotes that the 
ensemble was obtained by farming. A f denotes that the ensemble was farmed from a thermalised 
R algorithm dataset. We generate A'therm trajectories (of length r = 1/2) and a total of A^traj 
trajectories for each ensemble. 



Action (3 ^ Alg. St Mnemonic iV traj N therm 



DBW2 


0.72 


0.01/0.04 RHMC 


i 

50 


(D, 0.72, 0.01/0.04) 


6000 


1000 


DBW2 


0.72 


0.02/0.04 RHMC 


1 

54 


(D, 0.72, 0.02/0.04) 


6000 


1000 


DBW2 


0.72 


0.04/0.04 RHMC 


1 

50 


(D, 0.72, 0.04/0.04) 


3395 


1600 


DBW2* 


0.764 0.02/0.04 RHMC 


1 

70 


(D, 0.764, 0.02/0.04) 2940 


800 


DBW2*t 0.764 0.04/0.04 RHMC 


1 

70 


(D 0.764, 0.04/0.04) 


5320 


100 


DBW2 


0.78 


0.02/0.04 RHMC 


1 

70 


(D, 0.78 , 0.02/0.04) 


1505 


800 


DBW2 


0.78 


0.04/0.04 RHMC 


1 

70 


(D, 0.78 , 0.04/0.04) 


1620 


800 


Iwasaki 


2.13 


0.02/0.04 RHMC 


1 

50 


(I, 2.13 , 0.02/0.04) 


3595 


1000 


Iwasaki 


2.13 


0.04/0.04 RHMC 


1 

50 


(I, 2.13 , 0.04/0.04) 


3595 


1000 


Iwasaki 


2.2 


0.02/0.04 RHMC 


1 

50 


(I, 2.2 , 0.02/0.04) 


5900 


800 


Iwasaki* 


2.2 


0.04/0.04 RHMC 


1 

50 


(I, 2.2 , 0.04/0.04) 


5800 


800 



TABLE II: 2+1-flavor datasets used in our exploration of gauge actions with different plaquette 
and rectangle contributions and the ensemble mnemonics we use to describe them. The top three 
datasets have the same lattice spacing, and the others are finer. Here MD time refers to Molecular 
Dynamics evolution time. 

ci 3 ^ Alg. 5t Mnemonic MD time N therm 

-2.3 0.48 0.04/0.04 R 0.01 (C2.3, 0.48, 0.04/0.04, R) 800 300 

-3.57 0.32 0.04/0.04 R 0.01 (C3.57, 0.32, 0.04/0.04, R) 800 300 

-7.47 0.16 0.04/0.04 R 0.01 (C7.47, 0.16, 0.04/0.04, R) 760 300 

-3.57 0.333 0.04/0.04 R 0.01 (C3.57, 0.33, 0.04/0.04, R) 760 300 

-3.57 0.36 0.04/0.04 R 0.01 (C3.57, 0.36, 0.04/0.04, R) 690 300 

-2.3 0.53 0.04/0.04 R 0.01 (C2.3, 0.53, 0.04/0.04, R) 700 300 
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TABLE III: Parameter values for the 2+1 flavour RHMC algorithm simulations. A m i n and A max are 
the maximum and minimum eigenvalues of T>(rrii), which are needed for the rational approximation. 
nMD is the degree of the rational approximation used in the molecular dynamics evolution and umc 
is the rational approximation degree used in the Monte Carlo accept/reject step. The conjugate 
gradient stopping condition for the evolution was 1.0 x 10~ 6 and for the accept/reject step was 
1.0 x 10"" 10 . 

Ensemble mi m s mpy 





Amax 


Amin 




»MD n MC 


Amin 




»MD 


nMC 


Amin 




n MD "MO 


(D, 0.72, 0.01/0.04) 


2.4 


6 x 10" 


-5 


11 


18 


6 x 10" 


-l 


9 


14 


4 x 10" 


-2 


6 


9 


(D, 0.72, 0.02/0.04) 


2.4 


2 x 10" 


-4 


10 


15 


4 x 10" 


-l 


9 


14 


3 x 10" 


-2 


6 


9 


(D, 0.72, 0.04/0.04) 


2.4 


6 x 10" 


-4 


9 


14 


6 x 10" 


-l 


9 


14 


3 x 10" 


-2 


6 


9 


(D, 0.764, 0.02/0.04) 


2.42 


1 x 10" 


-4 


10 


15 


3 x 10" 


-i 


9 


14 


2 x 10" 


-2 


5 


8 


(D, 0.764, 0.04/0.04) 


2.42 


1 x 10" 


-4 


10 


15 


3 x 10" 


-4 


9 


14 


2 x 10" 


-2 


5 


8 


(D, 0.78, 0.02/0.04) 


2.42 


1 x 10" 


-4 


10 


15 


3 x 10" 


-1 


9 


14 


2 x 10" 


-2 


5 


8 


(D, 0.78, 0.04/0.04) 


2.42 


1 x 10" 


-4 


10 


15 


3 x 10" 


-1 


9 


14 


2 x 10" 


-2 


5 


8 


(I, 2.13, 0.02/0.04) 


2.4 


1 x 10" 


-4 


10 


15 


3 x 10" 


-1 


9 


14 


3 x 10" 


-2 


5 


8 


(I, 2.13, 0.04/0.04) 


2.4 


1 x 10" 


-1 


10 


15 


3 x 10" 


-1 


9 


14 


3 x 10" 


-2 


5 


8 


(I, 2.2, 0.02/0.04) 


2.4 


2 x 10" 


-4 


10 


15 


4 x 10" 


-1 


9 


14 


3 x 10" 


-2 


5 


8 


(I, 2.2, 0.04/0.04) 


2.4 


2 x 10" 


-4 


10 


15 


4 x 10" 


-1 


9 


14 


3 x 10" 


-2 


5 


8 



TABLE IV: Statistics for 2+1 flavour DBW2 ensembles with (3 = 0.72 generated with the RHMC 
algorithm. 

Ensemble < 5H > < eT m > Accept ~P^~ V (qq) (<?75<?) x 10~ 5 

(D, 0.72, 0.01/0.04) 0.388(14) 0.991(14) 0.658(6) 0.608201(21) 0.00315(1) 4.7(145) 
(D, 0.72, 0.02/0.04) 0.311(13) 1.000(15) 0.693(5) 0.608094(18) 0.00422(9) 9.6(7.9) 
(D, 0.72, 0.04/0.04) 0.398(16) 0.992(15) 0.651(9) 0.607788(10) 0.00612(11) -0.6(8.4) 
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TABLE V: Estimates of autocorrelation times for the average plaquette and the pseudoscalar 
meson on timeslice 12. Note plaquettes are separated by 1 trajectory while pseudoscalar meson 
correlators are measured on iV meas configurations each separated by 5 trajectories. 



Dataset Ensemble size 




u) C P (t = 


= 12) 




. m l \ AT . 
■ —J iV traj 


-^meas 


_cum 




T exp 


(D, 0.72 


, 0.01/0.04) 6000 


1000 


20(10) > 4 13(10) 


> 8 


(D, 0.72 


, 0.02/0.04) 6000 


1000 


10(5) 


> 5 15(10) > 12 


(I, 2.13 , 


0.02/0.04) 3595 


520 


6(3) 


> 6 5(3) 


> 5 


(I, 2.13 , 


0.04/0.04) 3595 


520 


6(4) 


> 4 7(3) 


> 9 
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TABLE VI: Results for the coefficients appearing in the fit to the static quark potential given in 
Eq. (|16p and the implied values of tq and 1/a. In computing the latter, the value tq = 0.5fm is 
used. The configurations used to obtain the results in this table are separated by 5 Monte Carlo 
steps. 



Ensemble Mnemonic 


■^v meas 


Vo a a 


ro 




1/a GeV 


(D, 0.72, 0.01/0.04) 


900 


0.880(10) 0.4bo(13) 0.0754(19) 


3.962(30) 


(16) 


1.564(12) 


(6) 


(D, 0.72, 0.02/0.04) 


o nn 

800 


0.853(11) 0.434(13) 0.0817(21) 


3.858(32) 


(75) 


1.523(13) 


(30) 


(D, 0.72, 0.04/0.04) 


280 


0.783(17) 0.365(22) 0.1018(32) 


3.554(31) 


(67) 


1.403(12) 


(27) 


(D, 0.764, 0.02/0.04) 


294 


0.816(20) 0.417(27) 0.0571(33) 


A n /I n / n 4 \ 

4.646(84) 


(21) 


1.833(33) 


(08) 


(D, 0.764, 0.04/0.04) 


300 


0.793(19) 0.398(26) 0.0648(31) 


4.397(60) 


(18) 


1.735(24) 


(07) 


(D, 0.78 , 0.02/0.04) 


180 


0.791(9) 0.412(11) 0.0524(16) 


4.863(54) 


(142) 1.929(21) 


(56) 


(D, 0.78 , 0.04/0.04) 


165 


0.761(10) 0.378(13) 0.0601(18) 


A r> r\r\ f A c\\ 

4.600(49) 


(16) 


1.815(19) 


(6) 


(I, 2.13 , 0.02/0.04) 


420 


0.834(10) 0.406(13) 0.0668(18) 


4.315(41) 


(97) 


1.703(16) 


(38) 


(I, 2.13 , 0.04/0.04) 


240 


0.792(14) 0.350(18) 0.0770(25) 


4.110(41) 


(21) 


1.622(16) 


(8) 


(I, 2.2 , 0.02/0.04) 


280 


0.804(7) 0.393(9) 0.0506(12) 


4.982(45) 


(30) 


1.966(18) 


( OQ\ 

[zo ) 


(I, 2.2 , 0.04/0.04) 


320 


0.781(6) 0.365(9) 0.0560(11) 


4.788(35) 


(30) 


1.890(14) 


(12) 


(I, 2.3 , 0.04/0.04) 


165 


0.741(5) 0.348(8) 0.0397(9) 


5.729(55) 


(137) 2.261(22) 


(54) 


(C2.3, 0.48 , 0.04/0.04, R) 


180 


0.814(16) 0.413(19) 0.0773(32) 


4.002(57) 


(60) 


1.579(22) 


(24) 


(C2.3, 0.53 , 0.04/0.04, R) 


160 


0.739(6) 0.384(8) 0.0492(12) 


5.074(48) 


(48) 


2.002(19) 


(19) 


(C3.57, 0.32 , 0.04/0.04, R) 


180 


0.828(15) 0.433(18) 0.0773(28) 


3.968(49) 


(43) 


1.566(19) 


(17) 


(C3.57, 0.33 , 0.04/0.04, R) 


145 


0.792(12) 0.416(15) 0.0657(21) 


4.335(47) 


(102) 1.711(19) 


(40) 


(C3.57, 0.36 , 0.04/0.04, R) 


135 


0.719(9) 0.371(12) 0.0500(17) 


5.060(66) 


(33) 


1.997(26) 


(13) 


(C7.47, 0.16 , 0.04/0.04, R) 


165 


0.832(15) 0.448(19) 0.0733(30) 


4.050(57) 


(430) 1.598(22) 


(17) 
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TABLE VII: Results for the lattice spacing extrapolated to the chiral limit for the two light quarks, 
m/ = — m res , while the mass of the third (strange) quark is held fixed at my = 0.04. The masses 
column lists those sea quark masses used in the extrapolation. 

Action (3 masses ro 1/a [GeV] 



DBW2 0.72 0.01, 0.02, 0.04 4.260(52)(12) 1.681(20)(5) 
DBW2 0.72 0.01, 0.02 4.177(115)(105) 1.648(41)(45) 
DBW2 0.764 0.02, 0.04 5.030(230) (040) 1.983(92)(16) 
DBW2 0.78 0.02, 0.04 5.184(134) (335) 2.046(53)(132) 
Iwasaki 2.13 0.02, 0.04 4.628(121) (239) 1.826(48) (94) 
Iwasaki 2.20 0.02, 0.04 5.239(114)(111) 2.068(44)(45) 



TABLE VIII: Results for pseudoscalar meson masses from fits to unitary RHMC data. 



Ensemble Mnemonic 


Fit range 






Results 










rap 


X 2 /dof 


correlators 


(D, 0.72 , 0.1/0.04) 


4-14 


0.01 


0.303(3) 


66/16 


SL-SL, SS-SS 


(D, 0.72 , 0.02/0.04) 


5-16 


0.02 


0.3742(9) 


34/18 


LL-LL, WL-WL 


(D, 0.72 , 0.04/0.04) 


5-16 


0.04 0.4916(10) 


50/18 


LL-LL, WL-WL 


(D, 0.764 , 0.02/0.04) 


7-15 


0.02 


0.311(1) 


20/12 


LL-LL, WL-WL 


(D, 0.764 , 0.04/0.04) 


6-16 


0.04 


0.4203(7) 


23/16 


LL-LL, WL-WL 


(D, 0.78 , 0.02/0.04) 


11-16 


0.02 


0.288(6) 


5/4 


WL-WL 


(D, 0.78 , 0.04/0.04) 


11-16 


0.04 


0.400(7) 


8/4 


LL-LL 


(I, 2.13 , 0.02/0.04) 


10-16 


0.02 


0.362(2) 


10/8 


LL-LL, SL-SL 


(I, 2.13 , 0.04/0.04) 


9-16 


0.04 


0.4665(9) 


14/10 


SL-SL j SS-SS 


(I, 2.2 , 0.02/0.04) 


5-16 


0.02 


0.315(2) 


29/18 


LL-LL, WL-WL 


(I, 2.2 , 0.04/0.04) 


6-16 


0.04 


0.425(1) 


21/16 


LL-LL, WL-WL 
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TABLE IX: Results for pseudo-scalar meson masses from fits to non-degenerate LL-LL DBW2 
(3 = 0.72 RHMC data. 





Fit range 


m { 




^min ^max 


0.01 0.02 0.04 


0.005 


8-16 


0.2619(26) 0.2769(30) 


0.01 


8-16 


0.3125(27) 


0.015 


8-16 


0.3321(23) 0.3450(25) 


0.02 


8-16 


0.3627(22) 0.3882(17) 


0.025 


8-16 


0.3914(21) 0.4031(22) 


0.03 


8-16 


0.4186(21) 0.4298(21) 0.4419(16) 


0.035 


8-16 


0.4445(20) 0.4552(20) 


0.04 


8-16 


0.4693(20) 0.4795(19) 



TABLE X: Results for the vector meson masses from fits to unitary RHMC data. 



Ensemble Mnemonic 


Fit range 






Results 




^min ^max 


my 


X 2 /dof 


correlators 


(D, 0.72 , 0.01/0.04) 


6-14 


0.580(10) 


13/12 


LL-LL, SL-SL 


(D, 0.72 , 0.02/0.04) 


6-16 


0.635(4) 


27/16 


LL-LL,WL-WL 


(D, 0.72 , 0.04/0.04) 


6-13 


0.703(5) 


11/10 


LL-LL, WL-WL 


(D, 0.764 , 0.02/0.04) 


8-14 


0.543(5) 


13/8 


LL-LL,WL-WL 


(D, 0.764 , 0.04/0.04) 


7-15 


0.607(3) 


21/12 


LL-LL,WL-WL 


(D, 0.78 , 0.02/0.04) 


10-15 


0.48(3) 


2/4 


WL-WL 


(D, 0.78 , 0.04/0.04) 


8-15 


0.575(5) 


24/13 


LL-LL,WL-WL 


(I, 2.13 , 0.02/0.04) 


5-14 


0.581(4) 


14/8 


SS-SS 


(I, 2.13 , 0.04/0.04) 


6-14 


0.661(3) 


21/12 


SS-SS, SL-SL 


(I, 2.2 , 0.02/0.04) 


7-16 


0.493(6) 


11/14 


LL-LL,WL-WL 


(I, 2.2 , 0.04/0.04) 


7-15 


0.586(5) 


19/12 


LL-LL,WL-WL 
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TABLE XI: Results for the vector meson masses from fits to non-degenerate LL-LL DBW2 (3 = 0.72 
RHMC data. 



1TL va \ 


Fit range 




mi 






^min ^max 


0.01 


0.02 


0.04 


0.005 


8-16 


0.574(22) 


0.586(18) 




0.01 


8-16 




0.596(15) 




0.015 


8-16 


0.595(13) 


0.610(13) 




0.02 


8-16 


0.610(11) 




0.648(12) 


0.025 


8-16 


0.6258(90) 


0.6405(97) 




0.03 


8-16 


0.6419(78) 


0.6564(87) 


0.6744(94) 


0.035 


8-16 


0.6579(69) 


0.6723(78) 




0.04 


8-16 


0.6740(62) 


0.6880(71) 





TABLE XII: Results for the nucleon and negative parity partner from fits to the unitary RHMC 
data. 



Dataset 




^min ^max 






Results 




Ensemble 


N from £Ibi 


N* from n Bl 


N* from Q.b 2 




m N * 


X 2 /dof 


(D, 0.72 , 0.02/0.04 ) 


10-16 


20-27 


6-11 


0.904(8) 


1.18(2) 


26/16 


(D, 0.72 , 0.04/0.04 ) 


8-14 


23-25 


4-8 


1.021(4) 


1.28(2) 


13/10 


(D, 0.764 , 0.02/0.04 ) 


11-15 


20-23 


6-10 


0.76(2) 


0.96(2) 


16/9 


(D, 0.764 , 0.04/0.04 ) 


10-16 


20-23 


7-12 


0.888(3) 


1.17(1) 


23/12 


( I, 2.13 , 0.02/0.04 ) 


9-14 


21-25 


7-9 


0.82(1) 


1.14(2) 


11/9 


( I, 2.13 , 0.04/0.04 ) 


8-12 


24-25 


6-11 


0.984(5) 


1.28(2) 


8/8 


( I, 2.2 , 0.02/0.04 ) 


8-16 


21-25 


6-10 


0.729(3) 


0.90(1) 


15/14 


( I, 2.2 , 0.04/0.04 ) 


10-15 


21-26 


6-11 


0.860(4) 1.051(5) 


19/13 
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TABLE XIII: Results for the residual mass for unitary data on the RHMC data. 



Ensemble Mnemonic Fit range Results 





^min ^max 


m[ es {m f ) 


X 2 /dof 


(D, 0.72 , 0.01/0.04) 


5-16 


0.01089(4) 


51/11 


(D, 0.72 , 0.02/0.04) 


9-16 


0.01092(7) 


18/15 


(D, 0.72 , 0.04/0.04) 


8-16 


0.01146(8) 


31/15 


(D, 0.764 , 0.02/0.04) 


9-15 


0.00535(2) 


18/13 


(D, 0.764 , 0.04/0.04) 


10-14 


0.00540(1) 


12/9 


(D, 0.78 , 0.02/0.04) 


11-15 


0.00428(2) 


6/4 


(D, 0.78 , 0.04/0.04) 


4-15 


0.00427(2) 


9/11 


(I, 2.13 , 0.02/0.04) 


10-15 


0.01127(3) 


12/11 


(I, 2.13 , 0.04/0.04) 


9-15 


0.01175(5) 


20/13 


(I, 2.2 , 0.02/0.04) 


10-15 


0.00688(2) 


19/11 


(I, 2.2 , 0.04/0.04) 


12-16 


0.00711(2) 


17/9 



TABLE XIV: Results for the residual mass, m' res (mf) for the LL non-degenerate DBW2 [3 = 0.72 
RHMC data. 



m val 


mf yn = 0.01 


mf yn = 0.02 


= 0.04 


0.005 


0.01101(6) 


0.01120(7) 




0.010 




0.01109(7) 




0.015 


0.01078(5) 


0.01097(6) 




0.020 


0.01068(5) 




0.01191(10) 


0.025 


0.01059(5) 


0.01077(6) 




0.030 


0.01050(5) 


0.01068(6) 


0.01167(9) 


0.035 


0.01042(5) 


0.01060(6) 




0.040 


0.01035(5) 


0.01053(6) 
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TABLE XV: The vector meson masses and the residual masses computed for the four data sets: 
(D, 0.72, 0.04/0.04, R), (C2.3, 0.48, 0.04/0.04, R), (C3.57, 0.32, 0.04/0.04, R), and (C7.47, 0.16, 
0.04/0.04, R). 

Ensemble m va i ray ^res( m /) 

(D, 0.72, 0.04/0.04, R) 0.02 0.626(10) 0.01141(6) 

0.03 0.655(8) 0.01119(6) 
0.04 0.685(7) 0.01100(5) 
0.05 0.715(6) 0.01084(5) 

(C2.3, 0.48, 0.04/0.04, R) 0.02 0.612(10) 0.00781(6) 

0.03 0.636(10) 0.00768(5) 
0.04 0.661(8) 0.00756(5) 
0.05 0.689(6) 0.00747(4) 

(C3.57, 0.32, 0.04/0.04, R) 0.02 0.598(11) 0.00754(5) 

0.03 0.628(9) 0.00742(5) 
0.04 0.657(8) 0.00732(5) 
0.05 0.686(7) 0.00723(4) 

(C7.47, 0.16, 0.04/0.04, R) 0.02 0.599(10) 0.00654(5) 

0.03 0.623(8) 0.00642(5) 
0.04 0.650(7) 0.00634(4) 
0.05 0.678(6) 0.00626(4) 
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TABLE XVI: Results for Z A from fits to unitary RHMC data. 



Ensemble Mnemonic 


Fit range 




Results 




((3,mi/m s ) 


^min £max 


Z A 


X 2 /dof correlators 


( D, 0.72 , 0.01/0.04 ) 


5-11 


0.7335(2) 


18/13 


SL-SL 


( D, 0.72 , 0.02/0.04 ) 


5-11 


0.7347(2) 


19/13 


SL-SL 


( D, 0.72 , 0.04/0.04 ) 


8-13 


0.7373(1) 


11/11 


SL-SL 


( D, 0.764 , 0.02/0.04 ) 


6-12 


0.75521(5) 


25/13 


WL-WL 


( D, 0.764 , 0.04/0.04 ) 


6-14 


0.75722(7) 


21/17 


WL-WL 


( D, 0.78 , 0.02/0.04 ) 


8-14 


0.7625(3) 


14/13 


WL-WL 


( D, 0.78 , 0.04/0.04 ) 


8-14 


0.7662(2) 


17/13 


WL-WL 


( I, 2.13 , 0.02/0.04 ) 


7-14 


0.73376(10) 


28/15 


SL-SL 


( I, 2.13 , 0.04/0.04 ) 


6-11 


0.7357(1) 


22/11 


SL-SL 


( I, 2.2 , 0.02/0.04 ) 


10-15 


0.74563(9) 


13/11 


WL-WL 


( I, 2.2 , 0.04/0.04 ) 


10-15 


0.74820(7) 


23/11 


WL-WL 



TABLE XVII: Results for Z A from fits to non-degenerate LL-LL DBW2 (3 = 0.72 RHMC data. 



™val 


Fit range 


mi 




^min ^max 


0.01 0.02 0.04 


0.005 


8-16 


0.73207(32) 0.73382(29) 


0.01 


8-16 


0.73413(26) 


0.015 


8-16 


0.73338(22) 0.73461(24) 


0.02 


8-16 


0.73416(20) 0.73518(29) 


0.025 


8-16 


0.73498(19) 0.73592(20) 


0.03 


8-16 


0.73583(19) 0.73668(20) 0.73651(22) 


0.035 


8-16 


0.73670(19) 0.73745(19) 


0.04 


8-16 


0.73759(19) 0.73825(18) 
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TABLE XVIII: Fitted values for the pseudoscalar meson decay constant using the value of Z4 and 
the axial-axial correlator for unitary RHMC data, Eq. (|40p . 



Dataset 


Fit range 




Results 






imin ^max 


Aa q a q 


fp 


X 2 /dof 


( D, 0.72 , 0.01/0.04 ) 


9-15 


0.028(2) 


0.099(3) 


9./6 


( D, 0.72 , 0.02/0.04 ) 


9-15 


0.044(1) 


0.111(2) 


11/6 


( D, 0.72 , 0.04/0.04 ) 


9-15 


0.075(3) 


0.128(3) 


11/6 


( D, 0.764 , 0.02/0.04 ) 


8-16 


0.0219(5) 


0.0890(8) 


14/7 


( D, 0.764 , 0.04/0.04 ) 

\ 1 if / 


8-15 


0.043(1) 


0.107(1) 


8/6 


( D, 0.78 , 0.02/0.04 ) 


11-16 


0.016(2) 


0.083(3) 


9/4 


( D, 0.78 , 0.04/0.04 ) 


12-16 


0.026(2) 


0.086(2) 


5/4 


( I, 2.13 , 0.02/0.04 ) 


5-16 


0.0307(7) 


0.096(1) 


13/11 


( I, 2.13 , 0.04/0.04 ) 


10-16 


0.057(2) 


0.114(2) 


3/6 


( I, 2.2 , 0.02/0.04 ) 


9-16 


0.0195(5) 0.0825(10) 


10/6 


( I, 2.2 , 0.04/0.04 ) 


10-16 


0.0407(7) 


0.1024(9) 


14/6 



TABLE XIX: Fitted values for the pseudoscalar meson decay constant using the value of Za and 
axial-axial correlator for non-degenerate DBW2 (3 = 0.72 RHMC data, Eq. (|40|) . 



m val 


mf yn = 0.01 


mf yn = 0.02 


mf yn = 0.04 


0.005 


0.0950(24) 


0.1017(23) 




0.010 




0.1043(22) 




0.015 


0.1009(22) 


0.1070(22) 




0.020 


0.1038(21) 




0.1213(29) 


0.025 


0.1066(21) 


0.1125(21) 




0.030 


0.1094(21) 


0.1152(21) 


0.1260(28) 


0.035 


0.1121(20) 


0.1180(21) 




0.040 


0.1147(20) 


0.1206(21) 
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TABLE XX: Fitted values for the pseudoscalar meson decay constant using the residual mass and 
the pseudoscalar density correlator, Eq. (|43p for the unitary RHMC data. 



Dataset 


Fit range 




Results 






^min £max 


Ap P 


fp 


X 2 /dof 


(D, 0.72 , 0.01/0.04 ) 


11-16 


0.76(5) 


0.097(3) 


18/5 


(D, 0.72 , 0.02/0.04 ) 


8-16 


0.93(3) 


0.115(2) 


7/8 


(D, 0.72 , 0.04/0.04 ) 


10-16 


0.96(4) 


0.129(3) 


9/5 


(D, 0.764 , 0.02/0.04 ] 


10-16 


0.49(1) 


0.0884(9) 


12/5 


(D, 0.764 , 0.04/0.04 ] 


11-16 


0.55(1) 


0.108(1) 


7/4 


(D, 0.78 , 0.02/0.04 ) 


11-16 


0.36(3) 


0.080(3) 


2/4 


(D, 0.78 0.04/0.04) 


11-16 


0.35(2) 


0.089(3) 


6/4 


(I, 2.13 , 0.02/0.04 ) 


10-16 


0.53(2) 


0.092(1) 


8/5 


(I, 2.13 , 0.04/0.04 ) 


9-16 


0.60(2) 


0.112(2) 


7/7 


(I, 2.2 , 0.02/0.04 ) 


9-16 


0.39(1) 


0.0833(8) 


10/7 


(I, 2.2 , 0.04/0.04 ) 


11-16 


0.48(1) 


0.1033(9) 


8/4 
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TABLE XXI: Fitted values for the pseudoscalar meson decay constant using Zj± and the local axial 
correlator, Eq. ([37]) . 



Dataset 


Fit range 


Results 






fp 


X 2 /dof 


( D, 0.72 , 0.01/0.04 ) 


10-16 


0.11(2) 


1/6 


( D, 0.72 , 0.02/0.04 ) 


10-16 


0.113(2) 


7/6 


( D, 0.72 , 0.04/0.04 ) 


12-16 


0.127(3) 


5/4 


( D, 0.764 , 0.02/0.04 ) 


11-16 


0.089(1) 


4/4 


( D, 0.764 , 0.04/0.04 ) 


11-16 


0.107(1) 


6/4 


( D, 0.78 , 0.02/0.04 ) 


12-16 


0.076(3) 


8/4 


( D, 0.78 , 0.04/0.04 ) 


10-16 


0.088(3) 


0.3/5 


( I, 2.13 , 0.02/0.04 ) 


11-16 


0.092(1) 


4/5 


( I, 2.13 , 0.04/0.04 ) 


10-16 


0.111(2) 


3/5 


( I, 2.2 , 0.02/0.04 ) 


11-16 


0.0830(8) 


7/5 


( I, 2.2 , 0.04/0.04 ) 


10-16 


0.102(1) 


10/5 
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TABLE XXII: Vector meson masses in the chiral limit and ro 



Dataset 


Results 


dataset 


my a" 1 (GeV) x 2 /do/ ro(fm) 



(D,0.72) 0.52±{ 1.49±| 1/1 0.564(13) 
(D.0.764) 0.45^ 1.71+J 0.580(20) 
(D.0.78) 0.38±| 2.1+2 _ q.487(57) 



(1,2.13) 0.46±} 1.68+4 0.544(34) 
(1,2.2) 0.37+] 2.08+J - 0.497(28) 



TABLE XXIII: Results for the residual mass and the vector meson masses with different rectangle 
contributions. The values shown are extrapolations to the valence chiral limit: m va i = for the 
residual masses and m va i = —m res for vector meson masses. The inverse lattice spacings determined 
from the vector meson masses are shown in the last column. 

Ensemble m ros m v V a (GeV) 



(D, 0.72, 0.04/0.04, R) 0.01176(7) 0.530(15) 1.450(41) 

(C2.3, 0.48, 0.04/0.04, R) 0.00802(6) 0.539(16) 1.430(43) 

(C3.57, 0.32, 0.04/0.04, R) 0.00773(6) 0.517(16) 1.489(45) 

(C7.47, 0.16, 0.04/0.04, R) 0.00670(6) 0.525(16) 1.467(44) 

(C3.57, 0.333, 0.04/0.04, R) 0.00480(9) 0.509(26) 1.51(8) 

(C3.57, 0.36, 0.04/0.04, R) 0.002306(21) 0.400(12) 1.92(7) 

(C2.3, 0.53, 0.04/0.04, R) 0.002683(30) 0.405(17) 1.90(8) 
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TABLE XXIV: Pseudo-scalar meson mass squared in the chiral limit and the bare s quark mass. 



Ensembles Results 

m P (m q -> 0) x 2 /dof m s (m p ) 

(D,0.72) -0.014(2) 0.006/1 0.049(2) 

(D.0.764) -0.004(2) - 0.044(3) 

(D.0.78) -0.011(10) - 0.036(5) 

(1,2.13) -0.002(4) - 0.042(2) 

(1,2.2 ) -0.008(3) - 0.032(2) 



TABLE XXV: Results for nix*, the nucleon and its negative partity partner. 

Dataset uik* "m-N rajy* 

(D.0.72) 0.588(7) 0.73(2) 1.02(5) 

(D.0.764) 0.51(1) 0.59(4) 0.80(3) 
(D.0.78) 0.46(2) 

(D.2.13) 0.529(10) 0.58(3) 0.92(5) 

(D,2.2) 0.443(9) 0.554(9) 0.71(3) 
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TABLE XXVI: Pseudoscalar meson decay constant in the chiral limit and corresponding lattice 
spacing from the physical value of /„-. 

Dataset Results 



(5 f PP a~ x (GeV) f£ A a~ l (GeV) ff a" 1 (GeV) 

0.72 0.080(5) 1.6(1) 0.082(4) 1.59(8) 0.082(4) 1.59(8) 

0.764 0.064(2) 2.04(6) 0.067(2) 1.94(7) 0.065(3) 1.99(9) 

0.78 0.068(8) 1.9(3) 0.07(1) 1.8(3) 0.062(8) 2.1(3) 

2.13 0.062(3) 2.10(9) 0.067(5) 1.9(1) 0.061(3) 2.1(1) 

2.2 0.057(2) 2.29(8) 0.056(3) 2.31(10) 0.057(3) 2.3(1) 
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P 

FIG. 1: Parameters of (3p and (3r for quenched (circles) and 2 flavor (squares) simulations with 
same lattice spacings and the choices for the parameters used for the simulations reported in this 
paper (diamonds). 
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t 

FIG. 2: Logarithm of the normalised autocorrelation function for the plaquette on the (D, 0.72 , 
0.02/0.04) dataset. 

i~exp is found from the slope at early t. 
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FIG. 3: Integrated autocorrelation time for the pseudoscalar meson on the DBW2 /3=0.72 datasets 
with the longest single Monte Carlo chains. The separation for decorrelated configurations should 
be 2r cum and the measurement is made every 5 trajectories. 
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FIG. 4: Values for vq found for four choices of the time variable t in Eq. (I14j) . These show that a 
reasonable plateau has been reached by t = 5, the value we adopt to determine r$. 
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(3=0.72 DBW2 RHMC 




FIG. 5: Chiral extrapolation of the Sommer parameter tq to the limit m Uj( j = mi + m rcs — > using 
results from the three DBW2, (3 = 0.72 ensembles with mi = 0.01,0.02 and 0.04. Extrapolations 
are shown using both the two lightest masses and all three. 
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FIG. 6: From top to bottom, the panels give the time history of the plaquette, (qq) and (q^/sq) and 
the toplogical charge for the (D, 0.72, 0.01/0.04) dataset. The values plotted are measured every 
5 time units. 



54 



0.0002 
0.0001 


-0.0001 
-0.0002 
-0.0003 
-0.0004 

0.0002 
le-04 


-0.0001 
-0.0002 

0.0002 
le-04 


-0.0001 
-0.0002 





0.0002 
0.0001 - 


-0.0001 



-0.0002. 







500 



1000 



1500 



2000 



2500 



3000 



FIG. 7: The time history of ((77517) is shown for the (D, 0.72, 0.01/0.04) dataset, with different 
sizes for the smoothing window. If a smoothing window of size s is used, the data plotted at time 
unit n is an average of data from n — s/2 to n + s/2 — 1. From top to bottom, the panels have a 
smoothing window of size 25, 50, 100 and 200 time units. 
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FIG. 8: The time history of the topological charge and (575(7), with a smoothing window of 50, is 
shown for the (D, 0.72, 0.01/0.04) dataset. The evolutions are very similar. 
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FIG. 9: Comparison of the 5Li and classically improved methods of calculating the topological 
charge for ~ 1000 HMC trajectory lengths on the (I, 2.13, 0.04/0.04) ensemble. 
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FIG. 10: The distribution of the topological charge for (from top to bottom) the (D, 0.72, 0.01/0. 04), 
(D, 0.72, 0.02/0. 04) and (D, 0.72, 0.04/0. 04) ensembles, taken from the classically improved method 
of calculating the topological charge. 
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FIG. 11: The evolution of topological charge for the four simulations (D, 0.72, 0.04/0.04, R) (top 
panel), (C2.3, 0.48, 0.04/0.04, R), (C3.57, 0.32, 0.04/0.04, R), and (C7.47, 0.16, 0.04/0.04, R) 
(bottom panel). 
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FIG. 12: The evolution of topological charge for the four simulations (C3.57, 0.32, 0.04/0.04, 
R)(top panel), (C3.57, 0.333, 0.04/0.04, R), (C3.57, 0.36, 0.04/0.04, R) and (C2.3, 0.53, 0.04/0.04, 
R) (bottom panel). Note that the topology stops evolving as we go to weaker couplings. 
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FIG. 13: Representative topological charge histories for the (D,0.764), (D,0.78), (1,2.13) and (1,2.2) 
actions (top to bottom). As can be seen, the rate of topological charge tunnelling decreases both 
when moving between the Iwasaki and DBW2 actions, and when moving to smaller lattice spacings. 



61 



0.5 

0.4 

0.3 

0.2 
0.4 

0.35 

0.3 

0.25 

0.2 
0.4 

0.35 
0.3 

0.25 
0.2 



1 I 1 I 1 I 1 



i i i ~r~ 



O LL-LL 

□ SL-LL 

O SL-SL 

A WL-WL 



J_ 



J_ 



A 



10 12 14 16 18 20 



i | i | i | i | i | i | i | i | i | r 



I 1 I 1 I 1 I 



J_ 



J_ 



A 



A 



$ m m $ 



J_ 



j_ 



a 



10 12 14 16 18 20 



A 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 
± A A ± A 



A 



A 



A 



A 



AAA 



A 



A 



A 



A A 



A 



A 



A 



2 4 6 8 10 12 14 16 18 20 
timeslice 



FIG. 14: The effective mass from the pseudoscalar meson correlator for the (D, 0.72, 0.01/0.04) 
ensemble with mj al = 0.01 for different sources. The top panel shows four different source/sink 
combinations, the middle panel is for a pseudoscalar meson correlator made of two SL quark 
propagators and the bottom panel is the WL-WL case. 
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FIG. 15: The vector meson effective mass for the (D, 0.72, 0.01/0.04) ensemble with mj al = 0.01 for 
different sources. The top panel shows four different source/sink combinations, the middle panel 
is for a vector meson correlator made of two SL quark propagators and the bottom panel is the 
WL-WL case. 
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FIG. 16: Comparison of smearing functions for a vector meson constructed from valence quarks 
with mass m = 0.04 using 72 Iwasaki (3 = 2.2 configurations with mi = 0.02 and m s = 0.04. 10 
configurations were averaged into each bin and then a full correlated analysis performed on the 
binned data. 
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FIG. 17: Nucleon effective mass plot for the (D, 0.72, 0.02/0.04) dataset. Circles correspond to 
the WL-WL-WL calculations, the squares for the SL-SL-SL calculations. 
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FIG. 18: The ratio R(t) that determines m res for the (D, 0.72, 0.01/0.04) ensemble, for pseudoscalars 
made from LL quark propagators and WL quark propagators. 



66 



0.755 



0.75 



0.745 



°- 74 



0.735 



0.73 



• LL-LL 

• SL-SL 



I 

if 



L 



a 725 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 l 

' 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 

timeslice 



FIG. 19: Za. for the (D, 0.72, 0.02/0.04) dataset. Different colors correspond to the different 
smearings. The lines shown are a fit to the SL-SL plateau. 
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FIG. 20: The ratio C A p{t)/C PP (t) versus time for the (I, 2.13, 0.02/0.04) and (I, 2.13, 0.04/0.04) 
datasets. The lines shown are the tanh fit to the LL-LL correlators. 
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FIG. 21: The m' res dependence on rrif for the DBW2 (3 = 0.72 dataset. The solid symbols show 
the unitary points used in linear extrapolation. 
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FIG. 24: f P from combined fits to NLO PQxPT for mp and f P on the (D, 0.72, 0.01/0.04) and 
(D, 0.72, 0.02/0.04) ensembles. The dashed symbols are excluded from the fits. 
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FIG. 25: The chiral extrapolation of the pseudoscalar meson mass squared versus quark mass for 
the two Iwasaki ensembles. The horizontal dashed lines show the physical kaon mass in lattice 
units as set by r$. 
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FIG. 26: The chiral extrapolation of the vector meson mass versus quark mass for the three DBW2 
ensembles. The vertical dashed lines show half the strange quark mass, enabling the m^ K (horizontal 
lines) to be predicted from each ensemble. 
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FIG. 27: The chiral extrapolation of the nucleoli (iV) and its negative parity partner (N*) for the 
DBW2 = 0.72 and = 0.764 datasets. 
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FIG. 28: The chiral extrapolation of the pseudoscalar decay constant versus quark mass for 
three methods. The upper panel shows the Iwasaki (3 = 2.13 dataset and the lower (3 = 2.2. 
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FIG. 29: Ratio of tuk* /mp versus (a/ro) 2 for all the datasets. The dotted lines are calculated from 
the ratio of the experimental values. 
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FIG. 31: The J parameter on all the datasets. 
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FIG. 32: Scaling of the baryon spectrum with lattice spacing squared. Closed circles denote the 
nucleon, N, and open circles the negative parity partner, N* at the chiral limit. Black symbols 
denote the experimental values scaled by appropriate factors of ro, red symbols the DBW2 ((3 = 
0.72,0.764) ensembles and blue symbols the Iwasaki {(3 = 2.13,2.2) ensembles. The value of 
ro = 0.5fm was chosen to give an indication of the experimental spectrum in these units. 
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FIG. 33: The Edinburgh plot. Red symbols denote the DBW2 ensembles and blue symbols the 
Iwasaki ensembles. The phenomenological curve derived from [551 ] has been shown to guide the 



eye. Experimental ratios and the values obtained in the static quark limit, where the hadron mass 
is equal to the sum of the valence quark masses, are given by the starred points. 
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FIG. 34: Ratio of fx/ 'fn for plotted against (a/ro) 2 for all the (3 values. The dotted lines are 
calculated from the ratio of the experimental values. 
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